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Introduction.- Description of granular rheology has 
been a long-term challenge for both science and technol¬ 
ogy. The problem extends to a vast range, from solid-like 
creep motion, gas-like, to liquid-like phenomena [1]. Sim¬ 
ilar to solid-liquid transitions, granular materials acquire 
rigidity when the density exceeds a critical value [2-5]. 
This phenomenon, referred to as the jamming transition, 
is universely observed in disordered materials such as col¬ 
loidal suspensions [6], emulsions, and foams [7], as well 
as granular materials. The jamming transition and its 
relation to the glass transition have attracted much in¬ 
terest in the last two decades, and various aspects have 
been revealed [8-12]. In particular, characteristics in the 
vicinity of the jamming point, including the critical scal¬ 
ing behavior, have been extensively investigated by ex¬ 
periments and numerical simulations [2-4, 13-25]. It has 
been shown that the shear stress, the pressure, and the 
granular temperature can be expressed by scaling func¬ 
tions with exponents near p ~ pj, where pj is the jam¬ 
ming transition density. The shear viscosity exhibits a 
form ~ (<pj — p)~ x with A ~ 2 for non-Browninan 
suspensions, foams, and emulsions [26-29], and a recent 
careful analysis demonstrated that A is in the range be¬ 
tween 1.67 and 2.55 [30]. It seems that the exponent A 
for granular flows takes a larger value than that for sus¬ 
pensions [18, 19, 25, 31], although a value in the range 
mentioned above has been reported as well [17]. How¬ 
ever, these studies are based on numerical simulations 
or phenomenologies without any foundation of a micro¬ 
scopic theory. 

Even when we focus only on the flow properties below 
the jamming point pj, which can be tracked back to Bag- 
nold’s work [32], we have not yet obtained a complete set 
in describing the rheological properties of dense granular 
flows. One of the remarkable achievements is the exten¬ 
sion of the Boltzmann-Enskog kinetic theory to inelastic 
hard disks and spheres [33-38]. However, it has been rec¬ 
ognized that the kinetic theory breaks down at densities 
with volume fraction p > pf = 0.49 [39-42], since there 
exists correlated motions of grains. A modification of the 


kinetic theory has been proposed [43], but a microscopic 
theory is still absent. 

Due to these situations, a microscopic liquid theory 
valid in the regime pp < p < pj has been desired. One 
attempt to respond to this problem is the extension of the 
mode-coupling theory (MCT) [44] for dense granular liq¬ 
uids. MCT has been applied to granular systems driven 
by Gaussian noises [45-47]. It qualitatively predicts a 
liquid-glass transition, though the noise in granular sys¬ 
tems is non-Gaussian in general [48-51]. MCT also has 
been applied to sheared dense granular systems [52-54]. 
There are three disadvantages of this approach: (i) the 
shift of p is necessary to describe the divergence of rj. 
(ii) Because of the shift of p, the jamming transition is 
not correlated with the divergence of the first peak of the 
radial distribution function, (iii) It predicts a plateau in 
the density correlation function, which is not observed in 
experiments nor in simulations [55-57]. 

From these observations, it is crucial to obtain an ex¬ 
plicit expression of the steady-state distribution function 
to construct a theory for dense granular liquids. For our 
purpose, we attempt to perform an expansion with re¬ 
spect to the dissipation to obtain an approximate explicit 
expression of the distribution function, valid in the weak 
dissipation regime for frictionless granular flows. Once 
the distribution is obtained explicitly, it is possible to 
calculate the steady-state average for arbitrary observ¬ 
ables. 

Microscopic starting equations.- We consider a three- 
dimensional system of N smooth granular particles of 
mass to and diameter d in a volume V subjected to sta¬ 
tionary shearing characterized by the shear-rate tensor 7. 
We assume that each granular particle is a soft-sphere, 
and the contact force acts only on the normal direction. 
For a simple uniform shear with velocity along the x axis 
and its gradient along the y axis, the shear-rate tensor 
is = A lS llx 5 uy = x,y,z) with a shear rate 7. It 
is postulated that the applied shear induces a homoge¬ 
neous streaming-velocity profile 7 r at position r, assum¬ 
ing that no heterogeneity such as shear banding exists. 
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Thus, the equation of motion is given by 

Pi =m(ri- 7 • n), (la) 

Pi = F^ +F[ vis) (lb) 

where F- cl) and Ff vls%> are, respectively, the elastic and 
the viscous contact forces acting on the grain i. Equa¬ 
tion (1) is known as the Shod equation, which is equiv¬ 
alent to Newton’s equation of motion under a uniform 
shear [58]. 

The most essential feature of granular systems, in con¬ 
trast to thermal systems, is that the steady state is de¬ 
termined by the balance between the viscous heating and 
the energy dissipation due to inelastic collisions. For 
sheared granular systems, this can be seen from the time 
derivative of the Hamiltonian, T-L(T) = YiLi P 2 /(2?n) + 
Yi j u(rij), where u(rij) is the inter-particle potential 
depending on = |r, — rj\, T = (r,;, p,;}Ei is the phase- 
space coordinate, and Yi j is the summation under the 
condition i ^ j. Then Ti(T) satisfies 

H(T) = -‘yV<rxv(T)-ZR,(T), (2) 

where 


cr^(r) 


i N 

f I" Pi,/J,Pi,v 

V ^ i m 


n,, (d?+dr’)] O) 


is the stress tensor and 

1 N 

n ^ = ~2 E^-^ (V1S) ( 4 ) 

Z i=1 


corresponds to the Rayleigh’s dissipation function [59, 
60]. For granular systems with the interparticle dis¬ 
sipative force proportional to the relative velocity, it 
is impossible to reduce the dynamics as overdamped. 


For later analysis, we assume that the contact forces 


deb 


= E 


and F™ - 


= E 


j^i ij 


F>- ’ are, re¬ 


spectively, given by F i 


(el) 


K,Q(d — Vij)(d — rij) and 


fY s) = -(0(d - nj)(vij ■ where 0 (x) = f 

for x > 0 and 0 (a:) = 0 otherwise, rij = r* — rj, 
fij = / tij , and Vij = Vi — Vj with the velocity of 


*-th particle t>,;. 

Steady-state distribution function.- To address the dis¬ 
tribution function for the nonequilibrium steady state, we 
start from an equilibrium state at t —> — oo and evolve 
the system with shear and dissipation. Then, the system 
is expected to reach a steady state at t = 0. Although it 
is impossible to derive an exact solution of the Liouvillc 
equation, equivalent to Eq. (1), for the 6 iV-dimensional 
distribution function, it is possible to obtain an approxi¬ 
mate solution by perturbation, parallel to the method for 
the linearized Boltzmann equation [61]. In the perturba¬ 
tion for dense sheared granular systems, it is simple to 
obtain the leading-order eigenfrequency of the relaxation 
towards the steady state [60]. Hence, we attempt to spec¬ 
ulate an approximate steady-state distribution, which we 


denote pgg(r), by applying an approximation which ex- 
plictly utilizes the relaxation time. 

For this purpose, we start from a formal but exact ex¬ 
pression for the distribution function [58], 


(ex) /y-A 

Pss ( r ) = ex P 


dr H eq (r(— t)) 


Peq(r(- 00 )), (5) 


which is the steady-state solution of the Liouvillc 
equation. Here, fi eq (r) = p e qH{T) — A(T) = 
—/3 eq [iVcrxyiT) + 2R,(T)\ — A(T) is the work function 
for pe q (r) = e -/3eq ^ r )/ f dr e~' 3e<1 ^ r ) at temperature 
Teq = where A(T) = {d/dT) • T is the phase-space 
volume contraction. We approximate Eq. (5) by intro¬ 
ducing the relaxation time T re i as 


exp 


dr H eq (r(— t)) 


p T r elfiss(r) 


( 6 ) 


which can be validated in the perturbation expansion 
of the Liouville equation around the canonical distribu¬ 
tion [60]. In the perturbation, we non-dimensionalize 
all the quantities, where the units of mass, length, and 
time are chosen as m, d, and \Jm/n, and introduce 
e = £/ y/Km < 1 as a perturbation parameter, which is 
related to the restitution coefficient e as e ~ y/2(l — e)/n 
for e ss 1. We attach a s tar * to the non-dimensionalized 
quantites, e.g. t* = t^Jn/m . Furthermore, we perform 
a scaling which leaves the steady-state temperature Tgg, 
which is the ensemble average of YiLiPi / fiNm) at the 
steady state in the dimensional unit, to be independent 
of e. This indicates that the granular fluid keeps its mo¬ 
tion in the limit e —> 0. From dimensional analysis, Tgg 
satisfies Tgs ~ ?n 3 d 2 7 4 /C 2 > which leads to T gg ~ e _ 2 7 * 4 . 
Hence, 7 * should satisfy 7 * ~ e 1 / 2 . We introduce a scaled 
shear rate 7 as 7 * = e 1 ^ 2 7 , where 7 is independent of e. 
We attach a tilde to the scaled quantities. The relax¬ 
ation time T re i is evaluated from the eigenfrequency of 
the perturbation expansion as 


Trei 



e we (Tgg) 


-1 


(7) 


in the hard-core limit [60], where we(T) = 
n^/T/m go(ip)d 2 is the Enskog frequency of 
collisions and go(<p) is the first-peak value of the radial 
distribution function. In Eq. ( 6 ), we have also introduced 


^ss(r) 


-Pss \iVv ( xy( T ) + 2A^ 1 g ) (r) 


( 8 ) 


where a 


(el) 

xy 


and Ai?gg* 


are respectively given by [60] 


Aftgg (T) = 7e ( 1 ) (T) + ^A(T), 


(9) 

(10) 

( 11 ) 
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Here, we ignore the contribution from the viscous shear 
stress, which is a higher-order correction in the limit e —> 
0. To summarize, we obtain 


g—-fss(r) 

PSs(r) = JdTe- 1 * s(r)’ 


( 12 ) 


where I SS (T) = ^ S H*{T) - f re iH ss (r) with H SS (r) = 
—/3 qs W*Vxy\r) + 2A^ 1 s ) (r) . We note that (i) the 

steady-state temperature Tss = Psi appears in Eq. (12), 

(ii) the steady-state average, (■ ■ ■ )ss = / drpss(r) • • • , 
is independent of the equilibrium temperature for p e q(r), 

(iii) the problem reduces to an ’’equilibrium” one with an 
effective Hamiltonian 7i e s(T) = Tss-fss(r) and the tem¬ 
perature Tss- Because the nonequilibrium contribution 
is small, we further expand pss(T) as 

e -« s w*(r) h+f rel f2 ss (r) 

pss(r) »- 1 (13) 


where gcs(p) = (1 — p/ 2 )/(l — p) 3 is the formula by 
Carnahan and Starling valid at p < pf [62], Note that 
we adopt the Kirkwood approximation for many-body 
correlations [60]. We further obtain the expression for 
the shear stress, 


<^j/( r )> 


ss ~ 


3y/6 ~2 S 3 / 2 

647r 3 / 2 ^ R 1 / 2 g 0 (p) 


(17) 


In the vicinity of the jamming point pj , S and R 
can be approximated as S « = 5^4 n* 3 go(p) 3 and R ss 
&' 3 n * 2 9o(p) 2 , respectively. This leads to the following 
expressions, 



(18) 

(19) 


from which we obtain 


with Z ~ f dT e Pss' H *( r ) 1 + r re iHss (r) . An approxi¬ 
mate expression for (A(T))ss is obtained as 

<A(r)) ss » <A(r)) eq + f rel (A(r)H ss (r)) e(i , (14) 

where (■■■)= f dT e~^ ssU T) ... i s the average with 
respect to the canonical distribution at Tss- It should 
be noted that Eq. (14) is the result of an exponential 
damping in the stress-stress correlation function in the 
Green-Kubo formula. 

So far Tgs is undetermined. We attempt to determine 
Tss by imposing the energy balance 


(^(r)) sg = -W (a xy (T)) ss - 2 (U(T)) SS = 0. (15) 


The explicit form of Tss will be given in Eq. (16). 

Shear viscosity and temperature.-Now we calculate the 
steady-state average of the shear stress and the energy 
dissipation rate by Eq. (14) and derive an explicit expres¬ 
sion for Tss- First, (cr a;y (r))ss is approximately given 

by <<7z„(r)> ss « -T re i7/?ss^* (T)^(T))^ . Sim- 


ilarly, the leading contribution gives 
(tZ^(T)) - 2f re i/3g S (^(^(rjA^s^r)) . Thus, we 

obtain the steady-state temperature from Eq. (D23) as 



where S and R are given by S = 1 + 

^ 2 n*g 0 (p) + y 3 n* 2 g 0 (ip) 2 + =5^4 n* 3 g 0 (p) 3 and 

R = n*go(ip) \38' 2 +n*go(p)\ , with P/i = 27t/15, 
2^3 = — 7 t 2 / 20 , ^4 = 3tt 3 /160, frt' 2 = -3/4, and 
= 77 t/16 [60]. We adopt the interpolation formula for 
hard spheres valid in the range pp<p<pj(pf = 0.49, 
pj = 0.639), g 0 (p) = gcs(Pf)(Pf - Pj)/{<P ~ <Pj), 


TSs-ipj-p)- 1 , ( 20 ) 

(da;y(r))gg ~ (pj — p) • (21) 

From Eq. (17), we obtain the shear viscosity g* = 
— (<M r ))ss/7 ~ (TJ - <P)~ 5/2 , or for rj' = 

-(^cnw^v^) cx -(^ y (r)) ss /i 2 , 

ij'~(pj-p)- 2 . (22) 

These results in Eq. (20)-(22) are consistent with the 
previous observations [26-30]. 

Comparison with simulation.- In order to verify the 
validity of the theory, we compare the theory with the 
molecular dynamics (MD) simulation. The parameters 
in the MD are N = 2000, e = 0.018375, and 7 * = 10“ 3 , 
10 -4 , 10 -5 . This condition corresponds to e = 0.96. 

The shear viscosity rj' and Tss are shown in Fig. 1, 
together with the results of the MD. We also show the 
log-log plots near pj as a function of pj — p in the inset. 
From the figure, we see that the theory agrees with the 
result of MD simulation for p < pj quantitatively with¬ 
out introducing any fitting parameter. The agreement is 
refined as the shear rate is decreased, where the hard-core 
limit is realized asymptotically. The smeared divergence 
in the vicinity of the jamming point observed for finite 
7 * is a well-known feature of the soft-core MD. We stress 
that the theory predicts the divergence of both the shear 
viscosity and the granular temperature as p —> pj, in 
contrast to the kinetic theory of inelastic spheres, where 
the shear viscosity behaves as rf ~ (pj — p)^ 1 and Tss 
remains finite [35]. On the other hand, in the MD, the 
shear viscosity behaves as rj' ~ (pj — p)~ 2 , in accordance 
with the theory. We note that the agreement between the 
theory and MD is relatively poor for Tss, though we have 
not clarified its reason. 

The relaxation time, Eq. (7), is shown in Fig. 2, to¬ 
gether with the result of the MD. The result of the MD is 
extracted from fitting by an exponential function for the 
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FIG. 1. (Color online) The density dependence of (a) the shear viscosity fj' and (b) the granular temperature. The result 
of the theory is shown in (blue) solid line, while that for the MD is shown in (red) triangles, (blue) diamonds, and (green) 
rectangles for 7 * = 10 — 3 ,10 — 4 ,10~ 5 . (Inset) The log-log plots for the results near ipj = 0.639. 



<P 


FIG. 2. (Color online) The density dependence of the re¬ 
laxation time, f re i = er* el . The result for the shear rate 
7 * = 1CP 4 is shown in (blue) solid line for the theory and 
(blue) diamonds for the MD simulation. 


transient data of the temperature relaxing to the steady- 
state. We see that Eq. (7) is quantitatively valid for 
ip < 0.63. 

Discussions.- From Eqs. (16) and (17), we see that the 
theory is subjected to the Bagnold scaling. The result of 
MD shows that the discrepancy from the Bagnold scal¬ 
ing becomes significant for ip > 0.635. Hence, there is 
room for improving the theory to cover the non-Bagnold 
regime. From the phenomenological scaling of jammed 
granules, the viscosity exhibits rj ~ |pj — ip\y<t>{ 1 - 2 /vi) ^ 
where y^ and y 1 are the scaling exponents for a xy ~ 
(ip - tpj) y * for ip > ip j and a xy ~ 7 ^ at ip = ipj [18]. 
If we assume y<$, = 1 as in Refs. [2, 3, 18, 19], we have 
y 1 = 4/7, which is close to the value of Ref. [17] [63] . For 
strongly dissipative situations, higher-order terms might 
alter the exponents of the divergences. Such a contribu¬ 
tion will be discussed elsewhere. 


Concluding Remarks. - We have developed a theory for 
jammed frictionless granular particles subjected to a uni¬ 
form shear with the aid of an approximate nonequilib¬ 
rium steady-state distribution function, and have shown 
that it remarkably agrees with the result of the MD simu¬ 
lation below the jamming point without introducing any 
fitting parameter. There are many future tasks for the 
application of our theory, such as the emergence of the 
shear modulus above the jamming point [2-5], the effect 
of friction for grains where the discontinuous shear thick¬ 
ening appears [64], the drag force acting on the pulling 
tracer [65-69], etc. Moreover, we should stress that the 
framework of our theory is quite generic. Indeed, we 
believe that the divergence of the viscosity for colloidal 
suspensions, 77 ~ (ipj — ip)~ 2 [70], can be understood by 
our framework. Therefore, the theory is expected to be 
applicable to a wide variety of phenomena in nonequilib¬ 
rium processes. 
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Otsuki and S. Takada for providing the prototype of the 
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supported by Grant-in-Aids for scientific research (Grant 
Nos. 25287098). The MD simulation for this work has 
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Institute for Theoretical Physics, Kyoto University. 


Appendix A: Microscopic starting equations 

We introduce the system we consider, i.e. a three- 
dimensional system of N smooth granular particles of 
mass to in a volume V subjected to stationary shearing 
characterized by the shear-rate tensor j Ml/ = 'yS MX 6 iyy . 
The time evolution of the system is determined by the 
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Newton’s equation of motion, 

mr t = i7 (el) + F j; (vis) (i = 1, • • • , iV), (Al) 

under a suitable boundary condition such as the Lees- 
Edwards boundary condition [58], accounting for the sta¬ 
tionary shearing. Here, r,; refers to the position of the 
i-th particle, and the dot denotes the time derivative. 
Fi is the conservative force, and is given 

by a sum Yl'j °f forces exerted on the *-th particle by 
other particles, 

F if ] = ' h ' n ' r '' 1 = ®(d - r ij)f(d - r^T-ij. (A2) 

Here, d is the diameter of the particle; u(r) is the pair 
potential; r tj = r* - rj, r tj = |r,y|, r i3 = r^/r^; and 
0(x) is the Heaviside’s step function, which is 1 for x > 0 
and 0 otherwise. Although a realistic elastic force might 
be Hertzian, where f{x) oc x 3 / 2 for a three-dimensional 
system, we adopt the linear spring model, f(x) = nx 
(k > 0 ), for simplicity. F^ V1S ' > = F-J 1 ^ denotes the 

viscous dissipative force which is represented by a sum of 
pairwise contact forces, 


rate of change of the internal energy H{T) can be calcu¬ 
lated, together with the Shod equation, Eq. (A4), as 


N 


i—1 

N 


« = £-■* + £' 


du(rij) 


ij 


dri 


Vi 


= E 


i=1 


P_ 

m l 
N 


F^+F^-i-Pi 


N 


= “7 


yE[— -+y-F, 

' L m 


(el)' 


-Ed"’ 

i= 1 
N 

• ^ m 1 


Pi 


l m 


+ 7 * ri 


i =1 
N 


Pi_ iji(vis) 
m 

N 




= -7 


= —}Va xy - 211. 


^(vis) 


(A 6 ) 


Here, a xy (F) is the cry-component of the stress tensor, 

1 N 

^(r) = yY. + (fg> + di -l )](AV) 


and 


F^ 1S> = -C ®(d - r ;i )r,, (■ Vij ■ r i3 )- (A3) 

Here, Vij = ly — Vj with Vi = ri refers to the rel¬ 
ative velocity of colliding particles, and ( (> 0 ) is a 
viscous constant corresponding to the harmonic poten¬ 
tial. This model for sheared granular systems has been 
widely adopted in computer-simulation studies. No¬ 
tice that the equation of motion is not invariant un¬ 
der the time-reversal map, since F^J 1S ^ changes sign for 
{ri,Vi} -a {r t ,-Vi}. 

Instead of Eq. (Al), we consider the following equation 
of motion valid for dense sheared granular systems, 

Pi = m (r» - 7 • n), (A4a) 

Pi = -F) (cl) + -F): (vIf>) - 7 • Pi, (A4b) 

which is called the Sllod equation [58]. Here, p, is the 
fluctuation of momentum around the steady-shear mo¬ 
mentum and referred to as the peculiar, or thermal, mo¬ 
mentum. The use of the Shod equation is restricted to 
uniform shear, which is in general not realized in realis¬ 
tic granular systems. However, it is a valid idealization 
in the higlr-density regime or for the flow on an inclined 
plane [40, 71, 72], where the profile of the shear velocity 
is well approximated as linear except for the region near 
the boundary. Hence, we adopt the Shod equation to 
address the rheology of dense sheared granular liquids. 

The internal energy of the system is given by 

«< r >=EfA+E'“< r«), (A5) 

2—1 i,j 

where T = {ri,Pi}f =1 is the phase-space coordinate and 
JA . is the summation under the condition i ^ j. The 


1 N i 

f cr) = - 2 E * • 7 V1S) = EE• 7T S) 

i=l i,j 

= i E ,0 ( d “ r ij)(vij • fij) 2 (A 8 ) 

corresponds to the Rayleigh’s dissipation function [59]. 
Note that we have utilized f/ v1s ^ = J2j^ i F,^ ls ' > and 

F}j ls) = —Fj} 1 ^ in the last equality. The first term on 
the right hand side of Eq. (A 6 ) is the work rate of shear 
exerted on the system and the second term is the energy 
dissipation rate due to inelastic collisions. In the steady 
state, the balance between these two are realized, i.e. 

H(T) = -}Va xy {T) - 21Z(T) = 0. (A9) 


Appendix B: Liouville equation 

The equation of motion, i.e. the Shod equation 
Eq. (A4), can be converted into the Liouville equa¬ 
tion [58]. For liard-sphere granular materials, it is con¬ 
ventional to adopt the pseudo-Liouville equation [73-75]. 
However, we start with the Liouville equation for soft 
spheres and take the hard-core limit later, since there is 
a subtlety in expressing the shear stress and the dissi¬ 
pation function via the pseudo-Liouvillian. We present 
the Liouville equation for uniformly sheared granular sys¬ 
tems in the following. It closely follows the formulation 
of Ref. [58], where similar explicit expressions can also 
be found in Refs. [76, 77]. 

The time evolution of phase variables, which do not 
explicitly depend on time and whose time dependence 
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comes solely from that of the phase-space point T = 
{ r i,Pi}iL i> is determined by 

j t Mnt)) = r • = iCA(T(t)). (B 1 ) 

The operator iC is referred to as the Liouvillian. For 
Eq. (A4), we have 


means of a perturbation expansion of Eq. (B7). In par¬ 
ticular, we attempt to construct a Rayleigh-Schrodinger 
perturbation theory, where the dissipation and the shear 
are treated as perturbations. 

In Eq. (B7), we consider the Laplace transform of 
P(r,i), i-e. 


iC = i£ (el) + iC.~ + i£ (vis) , 


(B2) 


,o 

’Mr) = / dte-^pi T,t), 

J — oo 


where the elastic part (i£( cl) ), the shear part (iC-y), and 
the viscous part ( iC ( V 1 S 1 ) are, respectively, given by 


which satisfies 


i£t(r)^„(r) = -z n ^ n {T). 


(Cl) 


(C2) 


«:"■> = E[- + F t" ~ 

L m ovi opi 

i 


Zj0"y - 


a - (7 • Pi) ' - 5 — 

or, dpi 


*£ (vis ) = y- F (vis) . 

^ * <9^ 


(B3) 

(B4) 

(B5) 


Here, n is an index for the Laplace modes, which is con¬ 
tinuous or discrete. This is an eigenvalue equation for the 
adjoint Liouvillian. The distribution function is given by 


p(r,t) = ^V”*Mr), 

n—0 


(C3) 


The formal solution to Eq. (Bl) can be written in terms 
of the propagator exp (iCt) as 


H(T(f)) = exp(i£t)A(r). 


(B 6 ) 


Hereafter, the absence of the argument t implies that 
associated quantities are evaluated at t = 0, e.g. T = 

W 

The Liouville equation for the nonequilibrium phase- 
space distribution function p(T, t) is given by 


where the summation over n is an integral for continuous 
modes. 

In order to perform a perturbation expansion, we first 
non-dimensionalize the observables by choosing the unit 
of mass, length, and time as to, d , and s/rn/n, and in¬ 
troduce an infinitesimal parameter 


e = 


c 


< 1 . 


'nm 


(C4) 


dp(T,t) 


d_ 

dt ~ ~ dT 

= —*£ f p(r, t), 


i>(M) 


r-£ + A, r) 


p(T,t) 

(B7) 


where *£t(T) = iC(T) + A(T) is referred to as the adjoint 
Liouvillian. Here, A(T) denotes the phase-space contrac¬ 
tion factor [58] which is defined by 

+ (B 8 ) 

For the Shod equation Eq. (A4), one finds 

A(r) = -— V' 0(d - nj) < 0. (B9) 

m 

The formal solution of the Liouville equation (B7) reads 
p(T, t ) = exp (-iCH) p(T, 0). (BIO) 


Appendix C: The perturbation expansion of the 
Liouville equation 


Note that the restitution coefficient e is related to e as 
e = exp [— £t c /m] via C,t c /m = n e/\J 2(1 — e 2 ), where 
t c = tt/ y/2n/ to — (C/m) 2 is the duration of contact of the 
spheres [31]. For e <C 1, the normalized energy dissipa¬ 
tion rate 1 — e 2 can be approximated as 1 — e 2 ss 2 £t c /m = 
V2ne + 0(e 3 ) ~ y/2ne. Together with 1 — e 2 ss 2(1 — e), 
we have e « \/2(l — e)/n for e ~ 1. We attach a star * 
to the non-dimensionalized quantites, e.g. t* = t\J k/to. 

Furthermore, we perform a scaling which leaves the 
steady-state granular temperature Tss, which is defined 
by / dr pss(r)]CfciP?/(3Abn), where pss(r) is the 
steady-state value of p(T, t), to be independent of e. This 
is because we are interested in the situation where the 
granular fluid keeps a meaningful fluid motion in the 
limit e —> 0. From dimensional analysis, (cr X!/ (r)}ss 
and (TZ{T)) SS scale as (o xy (T)) ss ~ jy/mTss/d 2 and 
(K(T)) SS ~ (Tss/m. From the energy balance, Eq. (A9), 
Tgs scales as Tgg ~ m 3 d 2 y 4 /C 2 , which leads to T gg ~ 
e~ 2 j* 4 . Thus, if we require T gg to be independent of e, 
7 * should scale as 


7 


A/ 2 


(C5) 


To be specific, we introduce a scaled shear rate 7 as 


• * 1/2 ~ 

7 = e ' 7, 


(C 6 ) 


We attempt to derive the eigenfrequencies of the dis¬ 
tribution function for dense sheared granular liquids by 


where 7 is independent of e. We attach a tilde to the 
scaled quantities. Although implicit, it should be kept in 
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mind that the anisotropic shear stress cr xy (T) is propor¬ 
tional to the shear rate. This implies that a xy (T) also 
scales as e 1 / 2 . To illuminate this feature, we introduce 
the following scaling for the anisotropic quantities, 

X T -'E xy -X = eV 2 X T • V xy • X, (C7) 

where X is an arbitrary anisotropic vector and 
C^xy)^ = >y For instance, we have P* tX P* iy = 
e 1/2 Pi,xPi,y 

Then, we can expand in terms of e as 


Explicitly, they are given by 

{ N N N 'i 

!, Ek*> E^ E^>> \ ■ ( C2 °) 

*=i »=i *=i j 

The equalities 

AT AT 

. £ (e q) *(r)^^ = £*£?* = 0 (Ai = *,!/, «)(C21) 

i=l i=l 

follow from the conservation of the momentum. We 
choose {0* (r)} (a = 1, • • • , 5) to be orthogonal, i.e. 


iC*\T) = *£( eq )*(r) + e*£i(r), (C8) 

where the unperturbed operator i£( eq )*(T) is given by 
Eq. (B3) and the perturbed operator iCi(r) reads 

i£i(T) =i£y(r)+i£ (vis) (r)+A(T) (C9) 

with 

N r g g - 

iA(r) = t£ , (cio) 

2=1 L ’ J 

N 

*£ (vis )(r) = E-^ (vis) • —, (cii) 

2—1 

^(vis) = _ £ e (1 _ r * )(r* • (C12) 

A(r) =-^'0( 1 -r* j ), (C13) 

i,3 

respectively. Accordingly, we expand the distribution 
function and the eigenvalue as 

KF) = p* eq ( t) [*<,°>*(r) + e^(r)] + o( e 2 ),(ci4) 

z* n = z£> + ezW+0(e 2 ), (C15) 

where p* q (r) i s the canonical distribution. Note that 
/9* q (r) satisfies 

j£ (eq )*(r)p* q (r) = o (cie) 

by virtue of the conservation of the internal energy. From 
Eqs. (C2), (C8), (C14), and (C15), we obtain 

*£ (eq )*(r)^)*(r) = _4°)*$(°)*(r), (ci7) 

,£( eq )* ( r)^a)(r) + p: q (r)- 1 i£ 1 (r) P : q (r)^°)*(r) 

= -4 0) *^ 1} (r) - 4 1 )^ 0) * ( r). (ci8) 

1. Eigenequation for the zero modes 

The unperturbed operator 'j£ (cq )*(r) has degener¬ 
ate five zero-modes, which we denote by </>*(T) (a = 

I,-" ,5), 

i£( eq )*(r)0*(r) = 0 (a = !,■■■, 5 ). (C19) 


J dr>: q (r)C(T)C'(r) cx <W- (C 22 ) 


Specifically, we adopt 

0 *(r) = i, 


1 f N r ,* 2 3 \ 

^a(r) = y- E \ - o NT * ’ ( C24 ) 

J%NT* \i= 1 1 1 J 


fi (r) = ipEi (C25) 

V 2=1 

where A = x, y, and z correspond to a = 3, 4, and 5, 
respectively. This leads to 

J dr* p e * q (r)C(r)C'(r) = 5aal . (C 26 ) 


As for the energy eigennrode, 4>2 (r), we consider only the 
kinetic energy. We expect this treatment to be valid for 
the unperturbed zero-modes in the hard-core limit. 

In the following, we work in the 5-dimensional space 
spanned by Eqs. (C23)-(C25). Then, the distribution 
function is explicitly given by 

P *(T,t) = j2^ + -pU(r) 

a=l 

x W*(T) + el* 1 ) (T) + • • • ] . (C27) 

Since {</>* (T)} (a = 1, • • • ,5) are five-fold degenerate, we 
must choose an appropriate linear combination to con¬ 
struct the unperturbed distribution function, 

*i 0 ) *(r) = E wC'(r), (C28) 

a'—l 

where {c af 3 }® /3=1 are the coefficients. 

We next consider Eq. (C18). Since Za'** = 0, we obtain 

= -4 1 V: q (r)4/W*(D. (C29) 

By multiplying Eq. (C29) by <j>* a , (T) and integrating with 
respect to T*, we obtain 

J dr*&,(r) [*£r(r)p* q (r)^w*(r) 

= -4 1 ) J dT* 4 >* a , (r) [pe q (r)^i 0) *(r)], (cso) 



which follows from the fact that the first term on the 
left-hand side of Eq. (C29) vanishes, 


J dT* <f>a>(T) [ P : q (r)*£^*(r)^ 1 )(r)] =o,(C3i) 

due to (r) = 0. From Eq. (C28), Eq. (C30) is 

expressed as 


^ zPc, 


^ \ Caot" W'a 


(C32) 


where 


w a . a „ = ydre(r)*£!(r)p: q (r)C4r). (C33) 


We obtain five eigenvalue equations for z„ ^ (a = 
1, • • • ,5) from Eq. (C32), 


Wc T = -z^c T 


(C34) 


where cj is a vector which constitutes a matrix c T , i.e. 
c T = {cf, ■■■ ,cj}. Here, the superscript T denotes a 
transpose. Then, the eigen equations read 


det 


W + z^l 


= 0 (a = !,•••, 5), (C35) 


where 1 is the identity matrix. 

Now we calculate the matrix elements of W, which can 
be decomposed into 


W = + VF (vis) + W (A) , (C36) 


where 


= 

aa' } 

[ dT* r a (T)t^(T) P : q (T)^(T), 

(C37) 

w (vi , s) = 

aa' f 

f dr* (r)i£ (vis) (r)p* q (r)^*,(r), 

(C38) 

W (A } = 

aa' ( 

[ dr*C(r)A(r) P : q (r)e(r). 

(C39) 


For Eqs. (C37) and (C38), it is necessary to evaluate the 
integrands, 

i£y(r)p e yr)C(r) 

= ; 4 (i> e * q (r)] C(r) + p: q (r)*£y( r )C(r), 

(C40) 

i 3 vis \T)p* eq (T)rjT) 

= bc (vis) (i>: q (r)] C(r) + p: q (r)^ (vis) (r)C(r). 

(C41) 


The Liouvillians act on p* q (T) as 

Pe q (r) _ 1 *A( r )Pe q (r) 

N 


d 


d 


Vi rPi,v a~ I 

OXi Opi'X J 


H*{T) 


= 

2=1 
N 

= FiY,[Pi'*Pi’V + Vi P iS ) \ =P*W*4f(T), (C42) 


2=1 


Pe q ( r ) _ 1 *£ (vis) (r)pe q (r) 

^ JL 

d P * 


N 


= -F E^ (vis) • 7^ n *W = -FY,Pi ■ ^ 


2=1 


2=1 


= \f E 0 ( 1 - r W (*« ■ *ii) (Pii ■ f ij) 

ij 

= \P* E 0 ( 1 - r tj ) ([P« + ~iVi e *\ ■ ru) ( P*j ■ fij) 


= 2/3*7t (1) (r) - e 1,2 i3*~iV*d ( Jy S){1) (X) 
« 2/3*7t (1) (T), 


(C43) 


where the scaled quantities axy (r), diy (T), and 
^W(r) are given by 


'xy 


i N 

iW> (r) = w £[: 


- . ~ pi(el) 

Pi,xPi,y + yiF>' 


2=1 


(C44) 


^ is)( 1 ) (r) = -^rE' (Pii • *v) 1 - r*j), (C45) 


*,3 


^ (1) ( r ) = zE'K-^) 20 (i-^)- 


On the other hand, the Liouvillians act on 0* (r) as 


*£y(r)C(r) = o (a = 1,4,5), 


(C46) 


N 


iCi(T)<%{ r) = - 


| NT* <= i 

N 


^ ^ Pi,xPi,y 


^(r)^(r) = -^=E«, 

and 

f£ (vis) (r)^(r) = 0 , 


*£ (vis) (r)^(r) = 


|N’T* i= i 


(C47) 
(C48) 

(C49) 

(C50) 


AT 


3 NT* 


(C51) 

AT p(vis) 

(a = 3,4,5), (C52) 

where A = x, y, and z correspond to a = 3, 4, and 5, 
respectively. Hence, Eqs. (C37) and (C38) are recasted 



9 


in the form 

W™, = J dr* p* q (r) 0 * (T)p**(V5-w (r)<j>* a , (r) 

+ J drv: q (r)C(r)*A(r)C^(r), (C53) 

w^ = J dr* Pe * q (r)C(r) 2 r^ ( 1 ) (r)C'(r) 

+ J dr>: q (r)C(r)*£^(r)^(r). (C54) 

Together with Eq. (C39), we obtain 

Waa' = J dT* p^(T)<t»l{T) [/3*W*vW(r) 

+2/3*A7?. (1) (r) +il^(T) +z£ (vis) (r)] T), (C55) 

where A7t (1 l(r) is given by 

A7^ (1) (r) = n w (r) + ^-A(r). (C56) 

We denote the three components of W as 

W = W (s) + W (al) + W^ 2 \ (C57) 

where 

w®, = J dr* P : q (r)C(r) 


x/3* 


^Vcr^ + 2A7t (1) l 0;, (T) (C58) 


is the symmetric part and 

w^) = J dr*p e * q (r)C(r)^(r)0*,(r), (C59) 

w< a J - / dr * Pe q ( r X(r)^ (vis) (r)</>^(r),(C60) 

are the asymmetric parts. The elements of VlA a b and 
Wt 3 - 2 ) are given from Eqs. (C47)-(C49) and Eqs. (C50)- 
(C52) as 


(C61) 



' 0 

0 

0 

0 

0- 


0 

0 

0 

0 

0 

FH( al ) = 

0 

0 

0 

0 

0 


0 

0 

-7 

0 

0 


. 0 

0 

0 

0 

0 . 


jy( a2 b 


0 0 0 0 


0 -~ c 

U 3 - 


0 0 0 
0 0 0 
0 0 0 
0 0 0 


(C62) 


where Sf is given by 

& = n* J d 3 r* g(r*,ip)Q{ 1 - r*). 


(C63) 


Here, g(r, <p) is the radial distribution function at equi¬ 
librium defined by ~ r i))e q = Nng(r,ip). In 

Eq. (C62), terms of order e 1 / 2 are neglected. The ele¬ 
ments of W-' s ) are given by 


jy( s ) = 


0 


0 

0 

0 


0 

/P* 

0 

0 

0 

ypf 


0 

0 

0 



& 

0 

0 

0 

0 

0 

2 ( $ xx 

j + 2 <s xy 

2 C S XZ 


0 

0 

If 

7 

0 

0 

0 

4 + 2fS yx 

2 yyy 

2 ( SV Z 


0 

0 

7 


0 

0 

0 

2<S ZX 

2 y zy 

2<S ZZ 


0 

0 

0 

0 



(C64) 


with = n* f d 3 r* g(r*,y)Q{ 1 - r*)r M f", where the 
anisotropic components vanish, = 0 (/1 / r), and 
the isotropic components are given by S ? xx = vv = ( S ZZ 
= Sf/3. From Eqs. (C57), (C61), (C62), and (C64), we 
obtain 


0 0 0 0 0 

^JlN& §Sf 0 0 0 

0 0 |Sf ^ 0 

0 0 0 §S? 0 

0 0 0 ‘0 §Sf 


(C65) 


Then, from Eq. (C35), the eigenvalues Za^ {a = 1, • • • ,5) 
are obtained as solutions of the following equation for A, 



0 
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clet 


fJVSf l^ + A 


f + X 7 0 

0 §Sf + A 0 
0 0 §S? + A 


= 0 . 


(C 66 ) 


(i) 


2. Eigenvalues of the first order z, 


From Eq. (C 66 ), the eigenvalues are given by 


;( 1 ) _ 


= o, 


z™ = 


(a = 2,3,4, 5), 


(C67) 

(C 68 ) 


from which we see that the four non-zero modes are de¬ 
generate and hence can be approximated as a single re¬ 
laxation mode. The relaxation time for this mode is given 

by 

n -1 

(C69) 


1 


'rel ' 


eSa ) 


In the hard-core limit, reduces to 

Sf-► V^4(T*), (C70) 

(cf. Eq. (Gil) which will be shown later), where uje(T) = 
Ai/tt n^/T/m go(ip)d 2 is the Enskog frequency of colli¬ 
sions, and go(<p) is the first-peak value of the radial dis¬ 
tribution function, i.e. < 70 (v 5 ) = g(d,<p). From Eqs. (C69) 
and (C70), we obtain 


Ael 


2y/l T 


1 -1 


eu* E (T*) 


(C71) 


n eq (r) in another form. Note that the shear-stress ten¬ 
sor a xy {T) decomposes into 


a xy (T) = aif(T) + ai™\T), 


(D3) 


where 


i N 

(r) ss 1 £ [ + Vi Fff} , (D4) 

V r -- ' L m ’ J 


i =1 
1 N 

J- X r -—(vis) 




(D5) 


From Ff is) = E' F J is) = - E' F^ (r) can be 


rewritten as 


*’< r > 


xy 


2V 


hJ 


— 2V ^ ' ( v ij ' rij) &ijyij r ij®{d r ij )■ (D6) 


From Vij = Pij/m + 7 • ry with 7 M „ = 7 S^S^y, cr^ Is) (r) 
can further be decomposed into 

4l is) (r) = 4: is) (1) (r) + v™ w (r), (D7) 


Appendix D: Steady-state distribution function 

We derive an approximate explicit expression for the 
nonequilibrium steady-state distribution function, with 
the aid of the relaxation time derived in Sec C. We start 
from an equilibrium state at t —> — 00 and evolve the 
system with shear and dissipation, reaching a nonequilib¬ 
rium steady state at t = 0. A formal but exact expression 
for the nonequilibrium steady-state distribution function 
is given by [58] 


(ex) /-—x 

Pss ( r ) = ex P 


[ dr f2 eq (r(—r)) 

J — OO 


which satisfies = 0. Here, 


Peq (r( OO)), 

(Dl) 


0 eq (r) = (3 eq H(r) - A(r) 

= —peq [ 7 Va xy {T) + 2K(T)} - A(r) (D 2 ) 

is the work function for the equilibrium distribution 
Peq(r) = e _/3eqW ( r )/ f dre _/3eq ^ r ) with the tempera¬ 
ture T eq = fief}- To proceed, it is convenient to cast 


where 


°xy S) ( 1 ) (r) = • hj) XijVijTijQid - ry), (D 8 ) 

ij 


4l iS) (2) ( r ) - E '444 0 ( d - rij). (D9) 


*0 


Similarly, 7£(r) defined in Eq. (A 8 ) can be decomposed 
into 


7Z(r) = 7Z W (r) + 7Z (2) (r) + rc (3) (r), (Dio) 


with 


7e ( 1 ) (T) = ^ • fy) 2 0(d- ry), (Dll) 

hj 


7Z (2) (r) = Y ■ rij) %ijVij r Y'ijQ(d ^ 2 j ) ? (D12 ) 

i,j 


^• (3) ( r ) = ^ ,S: ijyij r ij e ( d ~ r v ) • 


(D13) 


h3 
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One easily recognizes from these expressions the following 
equalities, 

TZ^{T) = ~jVa^ s ^ 1 \T), (D14) 

n^(T) = -lVaif^ 2 \r). (D15) 

Using these results, U eq (r) can be expressed as 


fi eq (r) = -/3 e , 


7 E a xy (r) - 7 ^< is)( 1 ) (r) + 2 An^(T) 

(D16) 


where we have defined 


A^W(r) = ^W(r) + ^A(r). 


(D17) 


We attempt to obtain an approximate expression of 
Eq. (Dl) in the form of the canonical distribution and its 
correction. The expression Eq. (Dl) is a product of the 
canonical term and the exponential of the time integral 
of the work function. From Appendix C, the exponential 
factor can be rewritten as 

r z-o i 


exp 


/ dr U eq (r(— r )) 


pTrel^SSlT) 


(D18) 


L J —oo 


where r re i is the relaxation time, Eq. (C71). In Eq. (D18), 
we have defined 


n ss (r) = -p ss 



Wvxy ls){1) (r) + 2 A^ 1 S ) (r) 


(D19) 

which is equivalent to fl eq (r) except that the inverse 
equilibrium temperature /? eq is replaced by its steady- 
state value /3ss via Eq. (A6). Accordingly, we replace 
/3 eq in the canonical term of Eq. (Dl) by /3ss- This gu- 
rantees the independence of the steady-state average 


where the scaled quantites are given by 


1 N 

^(E) = y7 E [pi,xPi,y + > 

i =1 

(D28) 

4y s)(1) ( r ) = 2V* E ■ r ^) x vVi3 r ij Q (l 

ij 

r ij )) (D29) 

^ (1) ( r ) = iE'H-^) 2 0 ( 

i,j 

(D30) 

A(r) = -^'0(l-rb). 

(D31) 




From this, the order of magnitude of the three terms in 
U ss (r) can be evaluated as 

.* yv (ei)* ( r) = e ~ W *dW(T), (D32) 

rV*a^*(T) = e 2 W*^ m {T), (D33) 

K {1) *{T) = eK {1 \T). (D34) 

We retain cr^*(r), ^^^(r) and neglect the higher- 

order term c 4 y Is ^E*(r), i.e. 


^ss(F) » -e^s[^d^(r) + 2A^ 1 s ) (r)l .(D35) 


As for the relaxation time, we introduce a scaled relax¬ 
ation time f re i with r* el = e _ 1 f re i, i.e. 


Trel 



(D36) 


<■ ■ ■ )ss = J dT pss(r)--- 

on /3 eq . As a result, we obtain 

e -is s(r) 

PSs(r) = /rire-UsT)’ 

where 


(D20) 


(D21) 


Tss(r) = /?ssH(r) - T rel u ss (r). (D22) 

The steady-state temperature T$s = is determined 
by the energy balance condition, which is given from 
Eq. (A9) as 


Note that f re i is related to 7 via Xg g , whose explicit ex¬ 
pression will be given later in Eq. (E 6 ). Then, Jss(r) is 
approximated as 

7 s * s (r) » (3* SS K*(T) - f re if 2 ss (r), (D37) 


where D ss (r) = -$ s lyU^CT) + 2AK^(T) 


i(i)/ 


. We 


treat f re iUss(r) as a correction, and expand pss(T) as 


pss(r) 


e -0ss«*(r) 


1 + T r elUss(r) 


z 


(D38) 


(n F)) gs = -W (a xy (T)) ss - 2 (K(T)) SS = 0. (D23) 

Next, we consider the scaling with respect to e and 
evaluate the order of magnitude of the terms in pss(T). 
The three terms in Eq. (D19) exhibit the following scaling 
properties, 

aif*(T) = e 1 / 2 dif(T), (D24) 

T^ is)(1) *( r ) = e 3 / 2 c r ^ is)( 1 ) (r), (D25) 

H {1) *(T) = eTZ {1 \T), (D26) 

A*(r)=eA(r), (D27) 


with Z 


J dT e _/ 3 ss' H *(r) 


1 + TVeiUss(r) 


. An approx¬ 


imate expression for the ensemble average of an observ¬ 
able A(r) with the weight pss(T) can be obtained from 
Ecp (D38) as 


<A(r)) ss » (A(r)) eq + f rel (A(r)u ss (r)) eq ,(D39) 

where (■■■)= f dT e _/ 3 ss' H T)... j s the average with 
respect to the canonical distribution at the temperature 
Tss- 
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Appendix E: Shear viscosity and temperature 

The steady-state average of the shear stress and 
the energy dissipation rate is obtained by the formula 
Eq. (D39), from which we can derive an explicit expres¬ 
sion for Tss- From the scaling arguments, the leading 
contribution to the shear stress comes from the elastic 
component a xy (T). Hence, (cr X!/ (r))ss is approximately 
given by 

(a xy (T)) ss « -f lel W ss V* . (El) 

This corresponds to the evaluation of the Green-Kubo 
formula by the relaxation time approximation in the cor¬ 
relation function. Similarly, the leading contribution to 
the energy dissipation rate comes from TZ^ (T) and hence 
we obtain 

(6(r)) gs » (^ (1) (r)) eq - 2 f re i^ s (^ (1) (r)A^ 1 s ) (r)) cq . 

(E2 ) eq 

In these expressions, the anisotropic terms (&iy (r)) eq 
and (7t ( 1 Hr)di® 1 ) (r)) eq are identically zero. We obtain 
the following results in the hard-core limit from straight¬ 
forward calculations, 


In the vicinity of the jamming point ipj, S and R 
can be approximated as 5 k n* 3 go(<p) 3 and R « 
M' 3 n* 2 go((p) 2 , respectively. This leads to the following 
approximate expressions, 




37 2 ^4 

1 

l 3 


^ -y, 97T 

ss ~ ^M n * 9o{tp) = 2240 


7 n*g 0 (ip), (E10) 


(d X j/( r )) S s « ~ 7Wqt\ 7 T ss /2n * 3 9o(<p) 2 


1280 
27t r 5 / 2 


10240V35 


7 V 7 / 2 5 o (<7 /2 . 


(Ell) 


In the course of the derivation, we have utilized 
Eq. (D36). Note that this expression agrees with the 
scaling from the dimensional analysis and obeys the Bag- 
nold scaling, (<7 xy (r))ss oc 7 . From these expressions, 
we obtain the scalings 


(E12) 

<**„(r)>ss ~ T* s 1 / 2 5o (^) 2 - (<fij - ^)“ 5/2 - (E13) 

From Eq. (E9), we obtain the expression for the shear 
viscosity 77* = -(a xy (T))ss/j, 


(£ (1) (r)) eq = ^NT£ s u%(T£ a ), 


(E3) V* 


27tt 5 / 2 
10240 v / 35 


7 n 


7/2 9o(v) 512 


(^j-^)- 5 / 2 , (E14) 


(r w (t)Atz^(t)\ - -nt£ U e ( Ts ?\ 
\ K ssy Veq 8 ss n*g 0 (<p) 

x \M 2 +&3n*9o{<p)\, 



v* 


n*T c 


*2 

SS 


or for if = -(d xy (r))ss/(7v / ^) a -(^*i/( r ))ss/7 2 , 


(E4) 


97T 2 

1280 ? 


*3 


go(v ) 2 ~ (<pj - <pY 


(E15) 


x [l+y 2 n*g 0 (ip ) + ^3 n* 2 g 0 (ip) 2 + ^4 n* 3 g 0 (ip) 3 ] , (E5) 

where M 2 = 1, M 3 = 37 r/ 4 , M 2 = 27t/15, < 5*3 = — 7 t 2 / 20 , 
and oS ^4 = 37r 3 /160. Explicit calculations of the equilib¬ 
rium correlations, Eqs. (E3)-(E5), are shown in Sec. F in 
detail. To be specific, Eq. (E3) is derived from Eq. (FI), 

Eq. (E4) from Eqs. (F10) and (Fll), and Eq. (E5) from 
(F14), (F22), and (F31). Here, the terms with coefficients 
Mi and ,M l are contributions of the i-body correlations, 
and the first term of Eq. (E5) is the contribution of the 
kinetic stress. From the energy balance Eq. (D23) of or¬ 
der e, i.e. — 7 F* (cr X y(r)) gs — 2 (iZ^T)^ = 0 , we obtain 

the steady-state temperature Tss as 

T* - S iFfil 

Tss "32 *R' (E6) 

where S and R are given by 

S = l + y 2 n*g 0 (ip) + y 3 n* 2 g 0 (ip) 2 + y A n* 3 g Q (ip) 3 , (E7) 

R = n*g 0 (tp ) [M’ 2 + ^3 n *9o{9 >)], (E8) 

with M' 2 = —3/4, M' 3 = 77t/ 16. We further obtain the 
expression for the shear stress, 




3V6 ~ 2 S 3 ' 2 

647t 3 / 2 ^ R 1 / 2 g 0 (ifi) 
(E9) 


Appendix F: Equilibrium correlations 

We derive the equilibrium correlations, Eqs. (E3)— 
(E5). In the course of the derivation, the hard-core limit 
is taken with the aid of Eqs. (G 6 ) and (G12). 

First, Eq. (E3) is calculated as 

(«" l(r) L = i(E'Wr < ‘«) 2 e( 1 -7)) 

\ i,3 / ea 



»^Vn*T s * s J d 3 r* g(r*)<d(l — r*) 

« ^JVT s *s^(T s * s ), (FI) 

where we have utilized ,■ j 'S(r* — r*j)) eq ~ Nn*g(r*) 
with g(r) the radial distribution function at equilibrium 
in the third equality and Eq. (G12) in the last equality. 
Next, we deal with Eq. (E4), 

(tZ {1) A= (tZ (1) TZ {1 ^ + ^-(tZ (1) A^ . (F2) 
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FIG. 3. (Color online) Schematic figure for the three-body 
correlation. The angle 9 r between r and r is in the range 
[rr/3, 7t] . 


The first term is decomposed as 

rp*2 

±G s_ / 

i,j l,k 

+8T s *| /X] 


(E'E' A K) A ( r «) 


eq 


+8T s *f ( X! '^ (1)Aii/ (r-* )^ (1) ^(r*.) ) , (F3) 


h3 


eq 


where we have defined (r*) = 0(1 — r*)r fi f 1 ’/ 4. 

Here, JT . ; is performed under the condition that any 
two pair of particles ( i,j,l ) is different. On the other 
hand, the second term on the right hand side of Eq. (F2) 
is given by 


where g( 3 \r, r') is the triple correlation function at equi¬ 
librium. The triple correlation is conventionally approx¬ 
imated by the Kirkwood’s superposition approximation 
as g( 3 \r,r') pa g(r)g(r')g(\r — r'|), where g(r) and g(r') 
are the radial correlations and g(\r — r'|) is the angular 
correlation [78] (cf. Fig. 3). In the present case, we are 
interested in the case where the spheres ( i , j) and (i, l ) 
are in contact. To ensure the connectivity, we insert a 
step function as 

5 (3) (r,r') ps g(r)Q(d - r)g(r’)Q(d - r')g(\r - r'|), (F6) 
which leads to 

/ eq 

Nr ,* 2 r°° r°° 

dr* J dr'* g{r*)g{r'*)r* 2 r '* 2 

X0(1 — r*) 2 0(l — r'*) 2 

x J dS J dS' (r ■ r’) 2 g(\r* — r'*\). (F7) 

Here, J dS and J dS' represent the angular integrations 
with respect to the solid angles of r and r', respectively. 
It should be noted that there is a subtlety in the integral 
of the angular correlation, g(\r* — r , *|), since spheres j 
and l are in contact when 9 r = n/3, where cos 0 r = r ■ f'. 
The radial distribution function g{r) = 1 + h(r), where 
h(r) is the total correlation function, consists of a 6- 
function (contact) contribution at r pa d and a regular 
contribution, which is approximately 1. Hence, it is rea¬ 
sonable to approximate h(r) by the 5-function contribu¬ 
tion, which is given by 


h(r) 


1 

AipS 



(F8) 


rp* 

i SS 

2 



rn*2 

^ss 

4 


E'E A(r*,)A(r; fc ) 



Hence, the four-body correlation cancels out and the two- 
and three-body correlations remain. 

The three-body correlation is calculated as 





x (f • f') 2 0(1 — r*)0(l — r'*) 

^y^r*/dVy 3 >( r-,C* )( r.ff 

x©(l — r*)0(l — r'*), 



(F5) 


with numerically fitted coefficients, A pa 3.43, B p a 1.45, 
and C pa 2.25 [79]. Here, 6 > 0 is defined by <p = ipj( 1 — 
5) 3 , which is approximated as S sa (ipj — for 

kp pa ipj. The angular integral is evaluated as (cf. Fig. 3) 


J dS J dS 1 (r • r'Y [1 + h(\r* — r* 


r 1/2 

= 8n 2 J d {cos 9 r ) cos 2 9. 

= 37r 2 + f dz (z + 1) 


1 + h 


^\/2(l — cos 9 r 

1 2 


l--(z+l) 2 h(z ), (F9) 


where h(z) = \6A/(z/5 + C) 4 + (z/6 + C) 2 ] /(4tp6). 

Note that 9 r is restricted to 9 r £ [7r/3,7r] due to the 
constraint that any two pair of particles is differ¬ 

ent. In approaching the jamming point, i.e. S — > 0, h{z) 
behaves as h(z) ~ [(5/^) 4 + (5/z) 2 ] —>■ 0 and hence 

does not contribute to the angular integral. The radial 
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FIG. 4. (Color online) Schematic figure for the four-body 
correlation. 


integrations can be carried out with the aid of Eq. (G12), 
from which we obtain 

\ij,l / eq 

(F10) 

The two-body correlation is calculated, similarly to 
Eq. (FI), as 

(J2'K {1) * iv (r*-) \ 

' / eq 

('Z'8{r*-r; j )\ ©(1 — r*) 2 


= — / dV 

16 ./ 


TVn* 

"TeT 


gM j w>* - (fid 


where Eq. (G12) is utilized in the last equality. This 
concludes the derivation of Eq. (E4). 

Finally, we deal with Eq. (E5). From the definition of 
the elastic stress tensor, Eq. (D28), we obtain 


/^(el)~(el)\ = —!— { jVTc 

\ ocy w xy j y *2 l £ 


*2 

SS 


1 


\ i,j l,k 


, (F12) 


eq. 


where the first term is the kinetic stress and the second 
term is the contact stress. The contact stress consists of 
three components, i.e. the two-, three-, and four-body 
correlations. 


First, the two-body correlation is calculated as 

l (E ” \ N 

\ id / eq 

= \nu* dr*g{r*)F( el) * 2 r* 4 J dSx 2 y 2 

c\ pOO 

= —Nn* / dr* g{r*)F^* 2 r* A 1 (F13) 

15 Jo 

where F ( el )*( r *) = -0(1 - r*)u*'{r*) is the magnitude 
of the elastic force with u*(r*) the pair potential, and 
f dS ■ ■ ■ is the integration with respect to the solid angle 
of f. From Eq. (G6), we we obtain the following expres¬ 
sion in the hard-core limit, 


1 ' „*2 P (eD- 2 \ „ ?l N T* 2 n*go(<p). (F14) 


-( Y'y* 2 F { ' 


1 


Next, the three-body correlation is calculated as 

1 / p( e, ) x * p( el ) x * 

2 { / j VijVik^ij x ik 
\i,j,k 


eq 


« l Nn * 2 f d 3 r*J d 3 r'*g( 3 \r*, r'*)r*r'* 
xxyx'y'F^*(r*)F^ l) *(r'*) 

1 pOO nOO 

w -Nn* 2 / dr* / dr'*g(r*)g(r'*)r* 3 r'* 3 

2 Jo Jo 

xF^*(r*)F^*{r'*)e(l - r*)0( 1 - r'*) 
x J dS J dS'xyx'y 1 g(\r* — r*'\), (F15) 

where Eq. (F6) is applied in the second equality. Simi¬ 
larly to Eq. (F7), the (5-function contribution of g(\r* — 
r*' |) vanishes in the angular integrals. Hence, we obtain 


J dS J dS' xyx'y' g(\r*-r*'\) 


/ 2 r 1 r 2 * r 2n 

d(cos9 r ) / d(cos8') / dc/> r / dcj>'xyx'y' 

-l J -1 Jo Jo 

(F16) 


7r 

10 ' 


Here, cos 9 r = r f' and 8', J)' are defined as 

x' = sin 9' cos <J>, (F17) 

y = sin 9' sin <f>' , (F18) 

z' = cos9'. (F19) 

This gives x and y in terms of 9 ’, <f>', 9 r , and <f> r as 
-/2 

- I sin 0„ cos — 

1 + z' 

(F20) 


x' “ \ x' y 

x = ( 1-I sin 9 r cos <b r -sin 9 r sin i 

\ 1 + z'J r 1 + z' 

+X,' COS 8 r , 

x'y' . n j. , (, y' 2 \ ■ n ■ m 

y = -— sin 9 r cos m r + 1-— sin 9 r sin <b r 

1 + z’ V 1 + 2 / 


+y' cos 9 r , 


(F21) 
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FIG. 5. (Color online) Schematic figure for the four-body correlation. 


where (f> r is the azimuthal angle of f in the spherical 
coordinate where r' points in the z-direction. Note that 
9 r is restricted to 9 r € [7t/3,7t], similarly to Eq. (F9). 
The radial integrations can be carried out with the aid 
of Eq. (G6), from which we obtain 


1 

2 



7r 

20 ' 


-^iVT s *fr* 2 


3o(^) 2 - 


(F22) 

Finally, the four-body correlation is calculated as 


k ' eq 

E '" * / * * \ r-r(el)x* 

VijWik - Vu) F ij F ik 


l 


/ eq 

~ ^Nn* 3 J d 3 r* j d 3 r'* j d 3 r"*g^(r*,r'*,r"*) 

xy*(y'* - y"*)x(x' - f")F( el )*(r*)F( el) *(|r'* - r"*\), 

(F23) 


where g^ (r ,r' ,r") is the quadruple correlation at equi¬ 
librium. Here, JT ■ ; k is performed under the con¬ 
dition that any two pair of particles ( i,j,l,k ) is differ¬ 
ent. We change the integration variable from r"* to 
r'"* = r"* — r'*, which leads to 



= i «»- 3 / <iV / dV* ybV"V 4 >(r*, r'* 

xy*y"'*xx"'F^*{r*)F^*{r"'*). 


r'"*) 

(F24) 


Similarly to Eq. (F6), we adopt the following approxima¬ 
tion, 


g( 4 )( r , r , ) r "') ~ g(r)Q(d — r)g(r')Q(d — r')g(r'")Q(d — r'") 
Xfl(|r - r'|) S (|r' + r'"|)s(|r' + r’" - r|), (F25) 


where the step function is inserted to ensure the con¬ 
nectivity of the spheres (i, j), (k, Z), and (*, k). Then, 
Eq. (F24) is approximated as 



xr * 3 r '* 2 r ///*3 p(el)* 

x0(l — r*)0(l — /*)©(! — r'"*) 


x dS dS' dS'"xyx"'y" 


xg(\r* - r'*\)g(\r'* + r”’*\)g(\r'* + r'"* - r*|). (F26) 


Here, f dS, f dS ', J dS'" represent the angular integra¬ 
tions with respect to the solid angles of f, r', r’", respec¬ 
tively. 

We first consider the radial integration in Eq. (F26). 
From Eq. (G6), integrations with r* and r'"* give 
Tgg go(ip) 2 . On the other hand, the integration with r'* 
gives goi^)- Hence, the radial integration is given by 

T^goiv) 3 . 

Next we consider the angular integration in Eq. (F26). 
Similarly to Eq. (F7), the 5-function contribution of 
the raidal correlations g(|r* — g(\r'* + and 

g{\r'* + r'"* — r*|) vanishes. Hence, it is given by 


j dS J dS' J dS'"xyx'"y" 

= J d(cos9') J d(cos9 r ) J d(cos9'") 


i-i 

r 2n 


pzix p p 

J d(j) J d(j) r J d(f>”' xyx"'y"', (F27) 


where cos 9 r = r ■ r', cos 9'" = r'" ■ r', and 0', <f>' are 
defined by Eqs. (F17)-(F19). Then, x and y are given 
by Eqs. (F20) and (F21), and x!" and y'" are given in 
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terms of Oft , ftft, 9', and ft by the replacement 0 r — > 9ft, 
ft ftft in Eqs. (F20) and (F21). For the case (a), 
where 9 r £ [7t/3,7t] and Oft £ [0,7t/ 3] (cf. Fig. 5), sphere 
j and l can move freely without interference in the ft and 
ftft direction, i.e. ft £ [0, 27r] and ftft £ [0, 27t]. Then, 
the angular integral is evaluated as 

/■■§ r 1 

J ftcosOr) J 1 dftosdft) f (Or, Oft) = - — ,(F28) 

where f(0 r ,0ft) = 47r 3 (3cos 2 0 r — l)(3cos 2 Oft — 1)/15 is 
the result of the integration with respect to O’ , ft, ft, 
and ftft. Similarly, for the case (b), where 0 r £ [27t/3,7t] 
and Oft £ [0,27t/ 3] (cf. Fig. 5), the angular integral is 
evaluated as 

r ~3 r 1 -^3 

J 1 d(cos0r) J id (co S 0ft)f(0 r ,0ft) = - —.(F29) 

As for the case (c), where 0 r and Oft vary in the range 
0 r £ [7t/3, 27t/3], Oft £ [7r/3 ,0 r ] (cf. Fig. 5), we obtain 

11 o 

f 2 f 2 Q-7T-3 

/ d(cos0r) d(cos0ft)f(0 r ,0ft) = —. (F30) 

J — \ J cos 6 r 

Hence, the angular integration is given by 37 t 3 /80. To¬ 
gether with the radial integration, we obtain 


1 

4 




' * * 7-1 (el)®* T-,(el)a;* 

VijVlkFij F ik 



3tt 3 

160 


iVTg|'"* 3 


n~g 0 (ft 3 . 
(F31) 


This concludes the derivation of Eq. (E5). 


Appendix G: Hard-core limit 

In systems of soft spheres with contact interactions, 
the Heaviside’s step function, 0(d — r), appears in the 
interparticle forces as in Eqs. (A2) and (A3). This step 
function reduces to a delta function, 8(d—r), in the hard¬ 
core limit. Here, we derive formulas valid in this limit. 

First, we derive a formula for the step function in 
the elastic force. It is convenient to introduce a func¬ 
tion which is referred to as the cavity distribution func¬ 
tion [78], 

y(r) = g(r) eft)- 1 , (Gl) 

where 

e ( r ) = e -/3ss«(r) ( G2 ) 

with u(r) = (k/ 2) <d(d — r)(d — r) 2 the pair potential of 
the elastic force. It should be noted that e(r) reduces 
to the step function in the hard-core limit, nd 2 /Xss —> 
oo, and hence its derivative e'(r) reduces to the delta 
function, e'(r) -A S(d — r), or e'(r*) <1(1 — r*) in the 

non-dinrensionalized form. Let us consider the hard-core 
limit of 

*g(r*)F^*r* 2 , (G3) 



where iA el )*(r*) = —0(1 — r*)u*'(r*) is the magnitude 
of the elastic force with u*(r*) the pair potential, for 
illustration. From u'(r)e(r) = — Tsse'(r), we obtain the 
following expression 


poo poo 

/ dr*g(r*)F^*r* 2 = -/ dr*g(r*)Q(l - r*)u*'(r*)r* 2 


J 0 


Jo 


poo 

= -T s * s / dr*g(r*)e'(r*)r* 2 ^ 

Jo 

This corresponds to the replacement 

F {el) (r)g(r) -> -T ss S(d - r)g 0 (ft), 


(G4) 

(G5) 


or, in non-dimensionalized form, 

F^*(r*)g(r*) -a -T s * s <5( 1 - r*)g 0 (ft. (G6) 


Next, we derive a formula for the step function in the 
viscous force Eq. (A3), which appears in e.g. Eqs. (B9), 
(D8), and (Dll). In accordance with the elastic force, it 
is convenient to introduce a variant of the cavity distri¬ 
bution function, 

2/1/2 M = 9(r) ei/ 2 (r) _1 , (G7) 

where 

e 1/2 (r) = e-V^ssvG). ( G8 ) 

Similarly to Eq. (G2), e\/ 2 (r) reduces to the step func¬ 
tion 0(d —r) in the hard-core limit, and hence its deriva¬ 
tive eC 2 (r) converges to the delta function S(d— r). We 
consider the following quantity, which has dimension of 
(time) -1 , for illustration, 

<A(r) >-q = —-••«)) 

~—eNn\l— [ d 3 r g(r)Q(d — r). (G9) 

V m J 

From / 2 (r) t/i/2(r) = ftftss K0(d — r)g(r), we obtain 


( A ( r ))eq = - eNn \l^\f^ J d3r 


= — eNn\ 



/Tss 


d re 1/2 (r)y 1/2 (r) 


-t -‘ineNnJ — gftftd 2 , 

V m 

which corresponds to the replacement 


(G10) 



- r) g(r) 


— S(d- r)g 0 (<p). 

m 


(GH) 


In the non-dimensionalized form, it is given by 


0(1 - r *)g(r*)^Z(r*)g 0 (<p), (G12) 


where 

HIO 5 VTSsa-r-) = 4 £™( V) w-n. (G13) 
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Appendix H: Molecular dynamics simulation 

We briefly describe the methods of the molecular dy¬ 
namics (MD) simulation we have performed. The unit 
of mass, length, and time is chosen as to, d, 1 /to/k, 
and we attach * to non-dimenzionalized quantities, e.g. 
t* = tsjn/m. The parameters of the simulation are the 
volume fraction <p, the shear rate 7*, the dissipation rate 
e = £/ y/nm, and the number of spheres N. The condi¬ 
tions are N = 2000, e = 0.018375, 7* = 10~ 3 ,10" 4 ,HT 5 , 
while (p is varied in the range 0.50 and 0.66. 

The governing equation is the Shod equation for uni¬ 
formly sheared systems, Eq. (A4). This equation is inte¬ 
grated numerically by the velocity Verlet algorithm. The 
time step of the calculation is chosen as At* = 0.01. We 
start with thermalizing the system at an initial temper¬ 
ature Ti n i in the absence of shear and dissipation, i.e. 
by setting 7 = 0 and _F( V1S ) = 0 in Eq. (A4). After 
thermalization, we switch on the shear and dissipation 
simultaneously, and evolve the system until the temper¬ 


ature T(t) = J2iLi P?:(t) 2 /(3-ZVm) is relaxed and starts to 
fluctuate around a steady value. We regard this behavior 
as a signal for reaching the steady state, and this steady 
value is identified as the steady-state temperature, T$s- 
Note that Xss is determined solely by the balance of en¬ 
ergy, i.e. 7 and £, and is independent of the choice of 
Tj n i. We extract the relaxation time r re i by fitting the re¬ 
laxation behavior of T(t) with an exponential function, 

g t/ Trel 

Then, the ensemble average of the shear stress around 
the steady state is measured by calculating 


1 N 

«( r )) ss = 


p(el)* p(vis)* 
^i.x ' ^i.x 


(HI) 

In order to suppress the statistical error, we sum the 
values of Eq. (HI), sampled at some interval in a single 
run, and divide the sum by the number of samples after 
the run is terminated. We have verified that the errors 
between independent runs are negligible. 
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